Observational bias in behavioural ecology: a case study and teaching tool

library(tidyverse) # tidy style coding
library(brms) # Bayesian models
library(loo) # for information criteria
library(tidybayes) # Bayesian aesthetics
library(MetBrewer) # colours
library(kableExtra) # tables
library(patchwork) # putting plots together
library(DT) # for search- and saveable tables
library(pander) # for simpler tables
library(png) # to load images
library(grid) # to plot images
library(ggdag) # to draw dags

Load in the data

We analyse data collected across three cohorts of students (2021-2023) enrolled in the third year School of Biosciences subject, Animal Behaviour, at The University of Melbourne. We hereafter refer to students as observers.

data <- read_csv("data/pigeon_data.csv") %>% 
  mutate(Student_ID = as.factor(Student_ID),
         Year = as.factor(Year),
         Foraging_prop = (Foraging_percentage / 100)) %>%
  filter(Year %in% c("2021", "2022", "2023"),
         Foraging_percentage != "NA",
         Primer_understood != "NA") %>%
  filter(Student_ID != "60" & Student_ID != "68" & Student_ID != "70" & Student_ID != "72") %>% # remove students that completed the task multiple times
  #Student_ID %in% c("60", "68", "70", "72")) %>% 
  select(-c(First_name, Surname, Peck_mean)) %>%  # remove names when ready
  left_join(
    read_csv("data/gender_data.csv") %>% 
      mutate(Student_ID = as.factor(Student_ID))
  ) %>% 
  rename(Observer_ID = Student_ID) %>% 
  select(Observer_ID, Year, Gender, everything())

data_peck <- 
  data %>%
  filter(Peck_rate_2 != "NA") %>% 
  pivot_longer(cols = Peck_rate_1:Peck_rate_2, names_to = "Trial",
               values_to = "Peck_rate")
  
# Create a function to build HTML searchable tables

my_data_table <- function(df){
  datatable(
    df, rownames=FALSE,
    autoHideNavigation = TRUE,
    extensions = c("Scroller",  "Buttons"),
    options = list(
      dom = 'Bfrtip',
      deferRender=TRUE,
      scrollX=TRUE, scrollY=400,
      scrollCollapse=TRUE,
      buttons =
        list('pageLength', 'colvis', 'csv', list(
          extend = 'pdf',
          pageSize = 'A4',
          orientation = 'landscape',
          filename = 'full_dataset')),
      pageLength = 78
    )
  )
}


my_data_table(data)

Column explanations

Observer_ID: unique, anonymised identifier for each observer.

Year: year that the experiment was conducted.

Gender: upon enrollment at The University of Melbourne, students are asked to indicate their title. We identified women as observers that answered “Miss” or “Ms” and Men as those who that answered “Mr”. Those with entirely missing entries were coded as “NA”.

Bias_treatment: the primer the observer received, where ‘satiated’ indicates that the observers were provided information prior to a trial that suggested pigeons were fell fed, whereas ‘hungry’ indicated that the pigeons were in poor condition and hungry.

Expectation: we asked the observers to indicate whether they thought the pigeons would be hungry or satiated. We included this question to test whether observers were appropriately primed by their bias treatment.

Primer_understood: did the observer’s expectation match the primer they received?

Foraging_percentage: the percentage of pigeons that observers estimated to be foraging over a 15 second period, while observing a large flock.

Peck_rate_1: the number of times a single chosen pigeon pecked the ground over a 15 second period.

Peck_rate_2: the number of times a second chosen pigeon pecked the ground over a 15 second period.

Foraging_prop: proportion of pigeons estimated to be foraging

\(~\)

Question 1: is estimation of foraging percentage affected by bias manipulations?

\(~\)

Exploring allocation versus expectation 1

\(~\)

We find that 22 of the 78 observers indicated a feeding motivation expectation opposite to that implied by the primer they were allocated.

This suggests that Expectation may be a better predictor of foraging estimation than allocated Bias_treatment. The relationship between these and all other variables that we expect to play a role in this system are depicted in Figure 1.

We explicitly assess the effect of treatment and expectation by fitting two models:

1. a model with allocated primer (Bias_treatment) as the predictor variable

2. a model with indicated hunger expectation (Expectation) as the predictor variable

\(~\)

gg_simple_dag <- function(d) {
  d %>% 
    ggplot(aes(x = x, y = y, xend = xend, yend = yend, colour = Variables)) +
    geom_dag_point() +
    scale_colour_manual(values = c("Included in model" = met.brewer("Hiroshige")[4], "Not included in model" = "grey80")) +
    geom_dag_text(color = met.brewer("Hiroshige")[7]) +
    geom_dag_edges() + 
    theme_dag()+
    theme(legend.position = "bottom",
          legend.title = element_blank())
}


observer_bias_dag <- dagify(EF ~ PEW + PEM + SB + TF,
        PEW ~ BT,
        PEM ~ BT,
        SB ~ PEW + PEM,
       labels = c("TF" = "True\n Foraging", 
                  "EF" = "Estimated\n Foraging",
                  "PEW" = "Prior\n Expectation Men",
                  "PEM" = "Prior\n Expectation Women",
                  "BT" = "Bias\n Treatment",
                  "SB" = "Selection\n Bias")) %>% 
  tidy_dagitty(seed = 5)
  
  
 observer_bias_dag <- left_join(observer_bias_dag$data, tibble(name = c("BT", "EF", "PEW", "PEM", "SB", "TF"),
                                   Variables = c("Included in model", "Included in model", "Included in model", "Included in model", "Not included in model", "Not included in model"))) %>% 
   
 gg_simple_dag()


observer_bias_dag

Figure 1. A direct acrylic diagram showing the flow of causation in our biological system. Our bias treatment (BT) was designed to affect the prior expectation (PE) of observer’s, a subset of which were women (PEW), while the remaining were men (PEM). Prior expectations may affect the estimated level of foraging (EF) directly, or more specifically through biased selection of particular foragers (SB). Estimated foraging is also affected by the true level of foraging (TF) carried out by the flock of pigeons in the footage. We hypothesised that a priori expectations of observer’s would affect their foraging estimates. The paths connecting coloured variables show that this can be tested by modelling the effect of bias treatment on estimated foraging, or alternatively by directly modelling the effect of prior expectation on estimated foraging.

\(~\)

Model 1: allocated primer

# First let's model the effect of bias treatment on foraging estimation 

foraging_model_treatment <- brm(Foraging_prop ~ 0 + Bias_treatment,
                                     data = data, family = Beta,
                                     prior = c(prior(normal(0, 1.5), class = b),
                                               prior(exponential(1), class = phi)),
                                     iter = 6000, warmup = 2000, chains = 4, cores = 4,
                                     control = list(adapt_delta = 0.8, max_treedepth = 10),
                                     seed = 1, file = "fits/foraging_model_treatment")

foraging_model_treatment <- add_criterion(foraging_model_treatment, criterion = "loo", file = "fits/foraging_model_treatment")

foraging_model_treatment

##  Family: beta 
##   Links: mu = logit; phi = identity 
## Formula: Foraging_prop ~ 0 + Bias_treatment 
##    Data: data (Number of observations: 78) 
##   Draws: 4 chains, each with iter = 6000; warmup = 2000; thin = 1;
##          total post-warmup draws = 16000
## 
## Population-Level Effects: 
##                        Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Bias_treatmentHungry      -0.47      0.14    -0.74    -0.20 1.00    14927
## Bias_treatmentSatiated    -0.47      0.15    -0.76    -0.17 1.00    14749
##                        Tail_ESS
## Bias_treatmentHungry      12186
## Bias_treatmentSatiated    11819
## 
## Family Specific Parameters: 
##     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi     4.21      0.61     3.11     5.48 1.00    13695    12202
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

foraging_model_treatment_gender <- brm(Foraging_prop ~ 0 + Gender * Bias_treatment,
                                     data = data, family = Beta,
                                     prior = c(prior(normal(0, 1.5), class = b),
                                               prior(exponential(1), class = phi)),
                                     iter = 6000, warmup = 2000, chains = 4, cores = 4,
                                     control = list(adapt_delta = 0.8, max_treedepth = 10),
                                     seed = 1, file = "fits/foraging_model_treatment_gender")

Table S1. Posterior estimates of the percentage of pigeon feeding rate, split by the primer observers were allocated.

new_data <- tibble(Bias_treatment = c("Hungry", "Satiated"))

new_data %>% 
  left_join(data %>% group_by(Bias_treatment) %>% summarise(`n observers` = n())) %>% 
  cbind(fitted(foraging_model_treatment, newdata = new_data, summary = T) %>% 
          as_tibble() %>% 
          mutate(across(1:4, ~ .x *100),
                 across(1:4, round, 2))) %>% 
  rename("Estimated proportion foraging" = Estimate,
         "Bias treatment" = Bias_treatment) %>% 
  pander(split.cell = 20, split.table = Inf)

Bias treatment	n observers	Estimated proportion foraging	Est.Error	Q2.5	Q97.5
Hungry	42	38.47	3.25	32.24	44.99
Satiated	36	38.59	3.51	31.85	45.65

or the gender analysis version

# add gender predictions

new_data_gender <- expand_grid(Bias_treatment = c("Hungry", "Satiated"),
                               Gender = c("Woman", "Man"))

new_data_gender %>% 
  left_join(data %>% group_by(Gender, Bias_treatment) %>% summarise(`n observers` = n()) %>% 
          ungroup() %>% filter(!is.na(Gender))) %>% 
  cbind(fitted(foraging_model_treatment_gender, newdata = new_data_gender, summary = T) %>% 
          as_tibble() %>% 
          mutate(across(1:4, ~ .x *100),
                 across(1:4, round, 2))) %>% 
  rename("Estimated proportion foraging" = Estimate,
         "Bias treatment" = Bias_treatment) %>% 
  pander(split.cell = 20, split.table = Inf)

Bias treatment	Gender	n observers	Estimated proportion foraging	Est.Error	Q2.5	Q97.5
Hungry	Woman	30	40.34	3.86	32.91	48.1
Hungry	Man	11	33.3	5.79	22.36	45.21
Satiated	Woman	27	38.99	4.03	31.26	47.1
Satiated	Man	9	36.94	6.68	24.27	50.4

Model 2: indicated expectation

# fit the same model, except using participant expectation rather than allocated bias treatment

foraging_model_expectation <- brm(Foraging_prop ~ 0 + Expectation,
                                     data = data, family = Beta,
                                     prior = c(prior(normal(0, 1.5), class = b),
                                               prior(exponential(1), class = phi)),
                                     iter = 6000, warmup = 2000, chains = 4, cores = 4,
                                     control = list(adapt_delta = 0.8, max_treedepth = 10),
                                     seed = 1, file = "fits/foraging_model_expectation")

foraging_model_expectation <- add_criterion(foraging_model_expectation, criterion = "loo", file = "fits/foraging_model_expectation")

foraging_model_expectation

##  Family: beta 
##   Links: mu = logit; phi = identity 
## Formula: Foraging_prop ~ 0 + Expectation 
##    Data: data (Number of observations: 78) 
##   Draws: 4 chains, each with iter = 6000; warmup = 2000; thin = 1;
##          total post-warmup draws = 16000
## 
## Population-Level Effects: 
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## ExpectationHungry      -0.40      0.13    -0.66    -0.14 1.00    15078    12139
## ExpectationSatiated    -0.56      0.15    -0.87    -0.26 1.00    15944    12233
## 
## Family Specific Parameters: 
##     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi     4.25      0.62     3.12     5.55 1.00    16348    12206
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

foraging_model_expectation_gender <- brm(Foraging_prop ~ 0 + Gender * Expectation,
                                     data = data, family = Beta,
                                     prior = c(prior(normal(0, 1.5), class = b),
                                               prior(exponential(1), class = phi)),
                                     iter = 6000, warmup = 2000, chains = 4, cores = 4,
                                     control = list(adapt_delta = 0.8, max_treedepth = 10),
                                     seed = 1, file = "fits/foraging_model_expectation_gender")

Table S2. Posterior estimates of the percentage of pigeons foraging, split by the actual expectation of observers.

new_data_2 <- tibble(Expectation = c("Hungry", "Satiated"))

new_data_2 %>% 
   left_join(data %>% group_by(Expectation) %>% summarise(`n observers` = n())) %>% 
  cbind(fitted(foraging_model_expectation, newdata = new_data_2, summary = T) %>% 
          as_tibble() %>% 
          mutate(across(1:4, ~ .x *100),
                 across(1:4, round, 2))) %>%  
  rename("Estimated proportion foraging" = Estimate,
         "Indicated expectation" = Expectation) %>% 
  pander(split.cell = 20, split.table = Inf)

Indicated expectation	n observers	Estimated proportion foraging	Est.Error	Q2.5	Q97.5
Hungry	44	40.2	3.19	34.08	46.62
Satiated	34	36.38	3.56	29.61	43.56

the gender analysis version

# add gender predictions

new_data_gender_2 <- expand_grid(Expectation = c("Hungry", "Satiated"),
                               Gender = c("Woman", "Man"))

new_data_gender_2 %>% 
  left_join(data %>% group_by(Expectation, Gender) %>% summarise(`n observers` = n()) %>% 
          ungroup() %>% filter(!is.na(Gender))) %>%
  cbind(fitted(foraging_model_expectation_gender, 
               newdata = new_data_gender_2, summary = T) %>% 
          as_tibble() %>% 
          mutate(across(1:4, ~ .x *100),
                 across(1:4, round, 2))) %>% 
  rename("Estimated proportion foraging" = Estimate,
         "Indicated expectation" = Expectation) %>% 
  pander(split.cell = 20, split.table = Inf)

Indicated expectation	Gender	n observers	Estimated proportion foraging	Est.Error	Q2.5	Q97.5
Hungry	Woman	32	42.7	3.75	35.42	50.11
Hungry	Man	12	34.41	5.58	23.91	45.81
Satiated	Woman	25	35.94	4.12	28.17	44.2
Satiated	Man	8	35.64	6.96	22.38	49.6

\(~\)

Build Figure 2a-d

\(~\)

Get posterior means and difference contrasts

# treatment model

draws_treatment <-
  as_draws_df(foraging_model_treatment) %>% 
  mutate(Hungry = inv_logit_scaled(b_Bias_treatmentHungry) *100,
         Satiated = inv_logit_scaled(b_Bias_treatmentSatiated)*100,
         diff_contrast = (Hungry - Satiated)) %>% 
  select(Hungry, Satiated, diff_contrast) %>% 
  pivot_longer(names_to = "Treatment", values_to = "Posterior_estimate", cols = 1:3) %>% 
  mutate(Predictor = "Allocated primer")

p2 <- 
  draws_treatment %>% 
  filter(Treatment != "diff_contrast") %>% 
  ggplot(aes(x = Treatment, y = Posterior_estimate)) + 
  stat_halfeye(aes(fill = Treatment), .width = c(0.66, 0.95), alpha = 0.9,
               point_interval = "median_qi", point_fill = "white",
               shape = 21, point_size = 4, stroke = 1.5) +
  scale_fill_manual(values = met.brewer("Hiroshige", 2)) +
  coord_flip(ylim = c(25, 55)) +
  #geom_vline(xintercept = 0, linetype = 2) +
  #scale_y_continuous(breaks = c(, 0, 1)) +
  xlab("Allocated primer") +
  ylab("Estimated % pigeons foraging") +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank(),
        text = element_text(size = 14))

p3 <-
  draws_treatment %>% 
  filter(Treatment == "diff_contrast") %>% 
  ggplot(aes(y = Posterior_estimate)) + 
   stat_halfeye(aes(fill = Treatment), .width = c(0.66, 0.95), alpha = 0.9,
               point_interval = "median_qi", point_fill = "white",
               shape = 21, point_size = 4, stroke = 1.5, scale =0.5) +
  scale_fill_manual(values = met.brewer("Hiroshige")[4]) +
  coord_flip(ylim = c(-20, 20)) +
  geom_hline(yintercept = 0, linetype = 2, linewidth = 0.75) +
  #scale_y_continuous(breaks = c(, 0, 1)) +
  xlab(NULL) +
  ylab("Hungry - Satiated difference\ncontrast (% points)") +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        text = element_text(size = 14))

# expectation model

draws_expectation <-
  as_draws_df(foraging_model_expectation) %>% 
  mutate(Hungry = inv_logit_scaled(b_ExpectationHungry) *100,
         Satiated = inv_logit_scaled(b_ExpectationSatiated)*100,
         diff_contrast = (Hungry - Satiated)) %>% 
  select(Hungry, Satiated, diff_contrast) %>% 
  pivot_longer(names_to = "Treatment", values_to = "Posterior_estimate", cols = 1:3) %>% 
  mutate(Predictor = "Indicated expectation")

p4 <-
  draws_expectation %>% 
  filter(Treatment != "diff_contrast") %>% 
  ggplot(aes(x = Treatment, y = Posterior_estimate)) + 
   stat_halfeye(aes(fill = Treatment), .width = c(0.66, 0.95), alpha = 0.9,
               point_interval = "median_qi", point_fill = "white",
               shape = 21, point_size = 4, stroke = 1.5) +
  scale_fill_manual(values = met.brewer("Hiroshige", 2)) +
  coord_flip(ylim = c(25, 55)) +
  #geom_vline(xintercept = 0, linetype = 2) +
  #scale_y_continuous(breaks = c(, 0, 1)) +
  xlab("Indicated expectation") +
  ylab("Estimated % pigeons foraging") +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank(),
        text = element_text(size = 14))

p5 <-
  draws_expectation %>% 
  filter(Treatment == "diff_contrast") %>% 
  ggplot(aes(y = Posterior_estimate)) + 
    stat_halfeye(aes(fill = Treatment), .width = c(0.66, 0.95), alpha = 0.9,
               point_interval = "median_qi", point_fill = "white",
               shape = 21, point_size = 4, stroke = 1.5, scale =0.5) +
  scale_fill_manual(values = met.brewer("Hiroshige")[4]) +
  coord_flip(ylim = c(-20, 20)) +
  geom_hline(yintercept = 0, linetype = 2, linewidth = 0.75) +
  #scale_y_continuous(breaks = c(, 0, 1)) +
  xlab(NULL) +
  ylab("Hungry - Satiated difference\ncontrast (% points)") +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        text = element_text(size = 14))

Now build the plots for the models that include gender

# treatment model

gender_treatment_draws <-
  fitted(foraging_model_treatment_gender, 
         newdata = new_data_gender, summary = F) %>% 
  as_tibble() %>% 
  rename(Hungry_Women = V1, Hungry_Men = V2, Satiated_Women = V3, Satiated_Men = V4) %>% 
  pivot_longer(cols = 1:4, names_to = "Group", values_to = "Posterior_estimate") %>% 
  separate(sep = "_", col = Group, into = c("Treatment", "Gender")) %>%
  mutate(Posterior_estimate = Posterior_estimate*100) 
  
calculate_all_the_diffs <-
  fitted(foraging_model_treatment_gender, 
         newdata = new_data_gender, summary = F) %>% 
  as_tibble() %>% 
  rename(Hungry_Woman = V1, Hungry_Man = V2, Satiated_Woman = V3, Satiated_Man = V4) %>% 
  mutate(Woman_h_s_diff = Hungry_Woman - Satiated_Woman,
         Man_h_s_diff = Hungry_Man - Satiated_Man,
         diff_diff = Woman_h_s_diff - Man_h_s_diff) %>% 
  select(contains("diff")) %>%
  rename(`H-S (women)` = Woman_h_s_diff,
         `H-S (men)` = Man_h_s_diff,
         `Interaction` = diff_diff) %>% 
  pivot_longer(cols = 1:3, names_to = "diff_contrast", values_to = "posterior_diff") %>% 
  mutate(posterior_diff = posterior_diff*100)
  

gp1 <- 
  gender_treatment_draws %>% 
  ggplot(aes(x = Gender, y = Posterior_estimate)) + 
    stat_slab(alpha = 0.8, shape = 21, aes(fill = Treatment)) +
  #stat_halfeye(.width = c(0.66, 0.95), alpha = 0.9,
   #            point_interval = "median_qi", point_fill = "white",
   #            shape = 21, point_size = 4, stroke = 1.5,
   #            fill = met.brewer("Hiroshige", 5)[2]) +
  #scale_fill_manual(values = met.brewer("Hiroshige", 2)) +
  scale_fill_manual(values = c(met.brewer("Hiroshige", 10)[4], met.brewer("Hiroshige", 10)[6])) +
  coord_flip()+#ylim = c(25, 55)) +
  #geom_vline(xintercept = 0, linetype = 2) +
  #scale_y_continuous(breaks = c(, 0, 1)) +
  labs(x = "Gender", y = "Estimated % pigeons foraging", fill = "Allocated\nprimer") +
  theme_bw() + 
  theme(legend.position = "bottom",
        panel.grid.minor = element_blank(),
        text = element_text(size = 12),
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 10))



gp2 <-
  calculate_all_the_diffs %>% 
  ggplot(aes(x = diff_contrast, y = posterior_diff)) + 
    #stat_slab(alpha = 0.9, shape = 21, aes(fill = Gender)) +
  stat_halfeye(.width = c(0.66, 0.95), alpha = 0.9,
               point_interval = "median_qi", point_fill = "white",
               shape = 21, point_size = 3, stroke = 1.5,
               fill = met.brewer("Hiroshige", 10)[5]) +
  #scale_fill_manual(values = met.brewer("Hiroshige", 2)) +
  coord_flip()+#ylim = c(25, 55)) +
   geom_hline(yintercept = 0, linetype = 2, linewidth = 0.75) +
  #geom_vline(xintercept = 0, linetype = 2) +
  #scale_y_continuous(breaks = c(, 0, 1)) +
  xlab("Difference contrast") +
  ylab("% points") +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank(),
        text = element_text(size = 12))

# expectation model

gender_expectation_draws <-
  fitted(foraging_model_expectation_gender, 
         newdata = new_data_gender_2, summary = F) %>% 
  as_tibble() %>% 
  rename(Hungry_Women = V1, Hungry_Men = V2, Satiated_Women = V3, Satiated_Men = V4) %>% 
  pivot_longer(cols = 1:4, names_to = "Group", values_to = "Posterior_estimate") %>% 
  separate(sep = "_", col = Group, into = c("Treatment", "Gender")) %>%
  mutate(Posterior_estimate = Posterior_estimate*100) 
  
calculate_all_the_diffs_2 <-
  fitted(foraging_model_expectation_gender, 
         newdata = new_data_gender_2, summary = F) %>% 
  as_tibble() %>% 
  rename(Hungry_Woman = V1, Hungry_Man = V2, Satiated_Woman = V3, Satiated_Man = V4) %>% 
  mutate(Woman_h_s_diff = Hungry_Woman - Satiated_Woman,
         Man_h_s_diff = Hungry_Man - Satiated_Man,
         diff_diff = Woman_h_s_diff - Man_h_s_diff) %>% 
  select(contains("diff")) %>%
  rename(`H-S (women)` = Woman_h_s_diff,
         `H-S (men)` = Man_h_s_diff,
         `Interaction` = diff_diff) %>% 
  pivot_longer(cols = 1:3, names_to = "diff_contrast", values_to = "posterior_diff") %>% 
  mutate(posterior_diff = posterior_diff*100)

gp3 <- 
  gender_expectation_draws %>% 
  ggplot(aes(x = Gender, y = Posterior_estimate)) + 
    stat_slab(alpha = 0.9, shape = 21, aes(fill = Treatment)) +
  #stat_halfeye(.width = c(0.66, 0.95), alpha = 0.9,
   #            point_interval = "median_qi", point_fill = "white",
   #            shape = 21, point_size = 4, stroke = 1.5,
   #            fill = met.brewer("Hiroshige", 5)[2]) +
  scale_fill_manual(values = met.brewer("Hiroshige", 2)) +
  coord_flip()+#ylim = c(25, 55)) +
  #geom_vline(xintercept = 0, linetype = 2) +
  #scale_y_continuous(breaks = c(, 0, 1)) +
  labs(x = "Gender", y = "Estimated % pigeons foraging", fill = "Indicated\nexpectation") +
  theme_bw() + 
    theme(legend.position = "bottom",
        panel.grid.minor = element_blank(),
        text = element_text(size = 12),
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 10))



gp4 <-
  calculate_all_the_diffs_2 %>% 
  ggplot(aes(x = diff_contrast, y = posterior_diff)) + 
    #stat_slab(alpha = 0.9, shape = 21, aes(fill = Gender)) +
  stat_halfeye(.width = c(0.66, 0.95), alpha = 0.9,
               point_interval = "median_qi", point_fill = "white",
               shape = 21, point_size = 3, stroke = 1.5,
               fill = met.brewer("Hiroshige", 10)[5]) +
  #scale_fill_manual(values = met.brewer("Hiroshige", 2)) +
  coord_flip()+#ylim = c(25, 55)) +
   geom_hline(yintercept = 0, linetype = 2, linewidth = 0.75) +
  #geom_vline(xintercept = 0, linetype = 2) +
  #scale_y_continuous(breaks = c(, 0, 1)) +
  xlab("Difference contrast") +
  ylab("% points") +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank(),
        text = element_text(size = 12))

\(~\)

Question 2: is estimation of feeding rate affected by bias manipulations?

\(~\)

We asked observers to count the number of pecks of the ground that a selected pigeon made over a 1 minute period. We use the number of pecks that occur as a measure of feeding rate.

We first estimated a baseline peck rate by observing 35 pigeons. To select the pigeons we observed, we split a still image of the foraging video (taken at time zero) into a 43 x 21 cell grid. From the 122 cells that contained pigeons, 40 were chosen by random number generation (see code chunk below). In the event that multiple pigeons were present in the cell, we selected the most prominent to observe. Five observations were discarded - three due to overlap of the same pigeon between cells that were selected by the random number generator, and two more as the pigeons left the field of view during the video and could no longer be tracked.

# curly brackets run all lines included within them

{set.seed(1) # so that sample produces a reproducible sequence 
  sample(1:122, 40, replace = FALSE)}

##  [1] 121  68  39   1  34  87  43  14  82  59  51  97  85  21 106  54  74   7  73
## [20]  79 110  37  89 101 118 100  44 103  33  84  35  70 108  42  38  20  28 117
## [39]  96  91

The selected pigeons for baseline observation are shown in the image below

img <- readPNG("pigeon_selection.png")
 grid.raster(img)

\(~\)

Estimating baseline feeding rate

\(~\)

Load in the data

baseline_data <- read_csv("data/baseline_peck_data.csv") %>% 
  select(1:5) %>%  # remove the comments column
  pivot_longer(cols = 4:5, names_to = "Observation", values_to = "Peck_rate") %>% 
  mutate(ID = as.factor(ID),
         Observation = str_remove(Observation, "Peck_count_")) %>% 
  rename(Pigeon_ID = ID) %>% 
  filter(!is.na(Peck_rate))

my_data_table <- function(df){
  datatable(
    df, rownames=FALSE,
    autoHideNavigation = TRUE,
    extensions = c("Scroller",  "Buttons"),
    options = list(
      dom = 'Bfrtip',
      deferRender=TRUE,
      scrollX=TRUE, scrollY=400,
      scrollCollapse=TRUE,
      buttons =
        list('pageLength', 'colvis', 'csv', list(
          extend = 'pdf',
          pageSize = 'A4',
          orientation = 'landscape',
          filename = 'baseline_dataset')),
      pageLength = 78
    )
  )
}


my_data_table(baseline_data)

• X and Y represent grid coordinates.

• Pigeon_ID identifies a specific pigeon

• Observation indicates whether this was the first or second scoring for a single pigeon. We scored each pigeon twice as distant pigeons were difficult to observe and to ensure that the correct pigeon was tracked throughout the minute of observation.

• Peck_rate is the number of times the ground was pecked over a minute of observation.

\(~\)

Fit a simple model to estimate median peck rate

baseline_peck_model_zi <-
  brm(Peck_rate ~ 1 + (1|Pigeon_ID),
      family = zero_inflated_negbinomial(), data = baseline_data,
      prior = c( prior(normal(0, 1.5), class = Intercept),
                 prior(exponential(1), class = sd),
                 prior(exponential(1), class = shape),
                 prior(exponential(1), class = zi)),
      chains = 4, cores = 4, warmup = 2000, iter = 6000,
      file = "fits/baseline_peck_model")

# wrangle the output

baseline_peck_predictions <-
  baseline_peck_model_zi %>% 
  as_draws_df() %>% 
  mutate(Baseline_estimate = exp(b_Intercept),
         peck_rate_sd = exp(sd_Pigeon_ID__Intercept)) %>% 
  select(Baseline_estimate, peck_rate_sd)

baseline_data %>% 
  distinct(Pigeon_ID) %>% 
  summarise(`n pigeons observed` = length(Pigeon_ID)) %>%
  bind_cols(
fitted(baseline_peck_model_zi, summary = T, re_formula = NA) %>% 
          as_tibble() %>%
          distinct(Estimate, .keep_all = T) %>% 
          mutate(across(1:4, round, 2)) %>% 
  rename(`Baseline median peck rate / per min` = Estimate)) %>%
  pander()

n pigeons observed	Baseline median peck rate / per min	Est.Error	Q2.5	Q97.5
35	1.32	0.51	0.54	2.53

\(~\)

Experimental data

\(~\)

Exploring allocation versus expectation 2.0

\(~\)

Once again, we expect that attention paid to and/or comprehension of the primer statement has a large effect on observers’ perception of pigeon foraging.

Lets again fit our two models:

1. a model with allocated primer (Bias_treatment) as the predictor variable

2. a model with indicated hunger expectation (Expectation) as the predictor variable

\(~\)

Model 1: allocated primer

# First let's model the effect of bias treatment on peck rate

peck_model_treatment <- 
  brm(Peck_rate ~ 0 + Bias_treatment + (1|Observer_ID),
      data = data_peck, family = negbinomial,
      prior = c(prior(normal(0, 1.5), class = b),
                prior(exponential(1), class = sd)),
      iter = 6000, warmup = 2000, chains = 4, cores = 4,
      control = list(adapt_delta = 0.9, max_treedepth = 12),
      seed = 1, file = "fits/peck_model_treatment")

peck_model_treatment <- add_criterion(peck_model_treatment, criterion = "loo", file = "fits/peck_model_treatment")
  
peck_model_treatment

##  Family: negbinomial 
##   Links: mu = log; shape = identity 
## Formula: Peck_rate ~ 0 + Bias_treatment + (1 | Observer_ID) 
##    Data: data_peck (Number of observations: 156) 
##   Draws: 4 chains, each with iter = 6000; warmup = 2000; thin = 1;
##          total post-warmup draws = 16000
## 
## Group-Level Effects: 
## ~Observer_ID (Number of levels: 78) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.21      0.15     0.01     0.54 1.00     4493     7497
## 
## Population-Level Effects: 
##                        Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Bias_treatmentHungry       2.32      0.13     2.06     2.57 1.00    14993
## Bias_treatmentSatiated     2.16      0.14     1.89     2.44 1.00    19527
##                        Tail_ESS
## Bias_treatmentHungry      10405
## Bias_treatmentSatiated    11645
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape     0.91      0.13     0.69     1.18 1.00    13608     8978
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

peck_model_treatment_gender <- 
  brm(Peck_rate ~ 0 + Gender * Bias_treatment + (1|Observer_ID),
      data = data_peck, family = negbinomial,
      prior = c(prior(normal(0, 1.5), class = b),
                prior(exponential(1), class = sd)),
      iter = 6000, warmup = 2000, chains = 4, cores = 4,
      control = list(adapt_delta = 0.9, max_treedepth = 12),
      seed = 1, file = "fits/peck_model_treatment_gender")

Table S3. The estimated peck rate of foraging pigeons, split by the primer observers were allocated.

new_data %>% 
  bind_cols(fitted(peck_model_treatment, newdata = new_data, summary = T, re_formula = NA) %>% 
              as_tibble() %>% 
              mutate(across(1:4, round, 2))) %>% 
  left_join(data_peck %>% group_by(Bias_treatment) %>% 
              distinct(Observer_ID) %>% summarise(`n pigeons observed` = n())) %>% 
  rename("Estimated peck rate" = Estimate,
         "Bias treatment" = Bias_treatment) %>% 
  bind_rows(fitted(baseline_peck_model_zi, summary = T, re_formula = NA) %>% 
              as_tibble() %>%
              distinct(Estimate, .keep_all = T) %>% 
              mutate(across(1:4, round, 2)) %>% 
              rename("Estimated peck rate" = Estimate) %>% 
              mutate(`Bias treatment` = "Baseline") %>% 
              bind_cols(baseline_data %>% 
                          distinct(Pigeon_ID) %>% 
                          summarise(`n pigeons observed` = length(Pigeon_ID)))) %>% 
  select(`Bias treatment`, `n pigeons observed`, everything()) %>% 
  pander(split.cell = 20, split.table = Inf)

Bias treatment	n pigeons observed	Estimated peck rate	Est.Error	Q2.5	Q97.5
Hungry	42	10.23	1.34	7.83	13.12
Satiated	36	8.74	1.24	6.59	11.44
Baseline	35	1.32	0.51	0.54	2.53

add the gender specific estimates

# add the gender table

new_data_gender %>% 
  bind_cols(fitted(peck_model_treatment_gender, 
               newdata = new_data_gender, summary = T, re_formula = NA) %>% 
          as_tibble() %>% 
          mutate(across(1:4, round, 2))) %>% 
  left_join(data_peck %>% group_by(Bias_treatment, Gender) %>% 
               distinct(Observer_ID) %>% summarise(`n pigeons observed` = n())) %>% 
    rename("Estimated peck rate" = Estimate,
         "Bias treatment" = Bias_treatment) %>% 
  bind_rows(fitted(baseline_peck_model_zi, summary = T, re_formula = NA) %>% 
              as_tibble() %>%
              distinct(Estimate, .keep_all = T) %>% 
              mutate(across(1:4, round, 2)) %>% 
              rename("Estimated peck rate" = Estimate) %>% 
              mutate(`Bias treatment` = "Baseline", Gender = "-") %>% 
              bind_cols(baseline_data %>% 
                          distinct(Pigeon_ID) %>% 
                          summarise(`n pigeons observed` = length(Pigeon_ID)))) %>% 
  select(`Bias treatment`, `n pigeons observed`, everything()) %>% 
  pander(split.cell = 20, split.table = Inf)

Bias treatment	n pigeons observed	Gender	Estimated peck rate	Est.Error	Q2.5	Q97.5
Hungry	30	Woman	10.84	1.68	7.91	14.51
Hungry	11	Man	8.52	2.25	4.95	13.61
Satiated	27	Woman	9.34	1.56	6.71	12.82
Satiated	9	Man	8.12	2.37	4.63	13.7
Baseline	35	-	1.32	0.51	0.54	2.53

Model 2: Indicated expectation

# fit the same model, except using participant expectation rather than allocated bias treatment

peck_model_expectation <- brm(Peck_rate ~ 0 + Expectation + (1|Observer_ID),
                                     data = data_peck, family = negbinomial,
                                     prior = c(prior(normal(0, 1.5), class = b),
                                               prior(exponential(1), class = sd)),
                                     iter = 6000, warmup = 2000, chains = 4, cores = 4,
                                     control = list(adapt_delta = 0.95, max_treedepth = 12),
                                     seed = 1, file = "fits/peck_model_expectation")

peck_model_expectation <- add_criterion(peck_model_expectation, criterion = "loo", file = "fits/peck_model_expectation")


peck_model_expectation

##  Family: negbinomial 
##   Links: mu = log; shape = identity 
## Formula: Peck_rate ~ 0 + Expectation + (1 | Observer_ID) 
##    Data: data_peck (Number of observations: 156) 
##   Draws: 4 chains, each with iter = 6000; warmup = 2000; thin = 1;
##          total post-warmup draws = 16000
## 
## Group-Level Effects: 
## ~Observer_ID (Number of levels: 78) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.18      0.13     0.01     0.48 1.00     4710     7228
## 
## Population-Level Effects: 
##                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## ExpectationHungry       2.41      0.12     2.18     2.66 1.00    16705    11820
## ExpectationSatiated     2.00      0.14     1.73     2.28 1.00    18644    12176
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape     0.93      0.12     0.71     1.20 1.00    14884    11083
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

#loo_compare(peck_model_treatment, peck_model_expectation)

peck_model_expectation_gender <- brm(Peck_rate ~ 0 + Gender * Expectation + (1|Observer_ID),
                                     data = data_peck, family = negbinomial,
                                     prior = c(prior(normal(0, 1.5), class = b),
                                               prior(exponential(1), class = sd)),
                                     iter = 6000, warmup = 2000, chains = 4, cores = 4,
                                     control = list(adapt_delta = 0.95, max_treedepth = 12),
                                     seed = 1, file = "fits/peck_model_expectation_gender")

Table S4. The estimated peck rate of foraging pigeons, split by the indicated expectation of the observers.

new_data_2 %>% 
  bind_cols(fitted(peck_model_expectation, newdata = new_data_2, summary = T, re_formula = NA) %>% 
              as_tibble() %>% 
              mutate(across(1:4, round, 2))) %>% 
  left_join(data_peck %>% group_by(Expectation) %>% 
              distinct(Observer_ID) %>% summarise(`n pigeons observed` = n())) %>% 
  rename("Estimated peck rate" = Estimate,
         "Indicated expectation" = Expectation) %>% 
  bind_rows(fitted(baseline_peck_model_zi, summary = T, re_formula = NA) %>% 
              as_tibble() %>%
              distinct(Estimate, .keep_all = T) %>% 
              mutate(across(1:4, round, 2)) %>% 
              rename("Estimated peck rate" = Estimate) %>% 
              mutate(`Indicated expectation` = "Baseline") %>% 
              bind_cols(baseline_data %>% 
                          distinct(Pigeon_ID) %>% 
                          summarise(`n pigeons observed` = length(Pigeon_ID)))) %>% 
  select(`Indicated expectation`, `n pigeons observed`, everything()) %>% 
  pander(split.cell = 20, split.table = Inf)

Indicated expectation	n pigeons observed	Estimated peck rate	Est.Error	Q2.5	Q97.5
Hungry	44	11.27	1.4	8.81	14.29
Satiated	34	7.46	1.06	5.62	9.81
Baseline	35	1.32	0.51	0.54	2.53

add the gender specific estimates

# add the gender table

new_data_gender_2 %>% 
  bind_cols(fitted(peck_model_expectation_gender, 
               newdata = new_data_gender_2, summary = T, re_formula = NA) %>% 
          as_tibble() %>% 
          mutate(across(1:4, round, 2))) %>% 
  left_join(data_peck %>% group_by(Expectation, Gender) %>% 
               distinct(Observer_ID) %>% summarise(`n pigeons observed` = n())) %>% 
    rename("Estimated peck rate" = Estimate,
         "Indicated expectation" = Expectation) %>% 
  bind_rows(fitted(baseline_peck_model_zi, summary = T, re_formula = NA) %>% 
              as_tibble() %>%
              distinct(Estimate, .keep_all = T) %>% 
              mutate(across(1:4, round, 2)) %>% 
              rename("Estimated peck rate" = Estimate) %>% 
              mutate(`Indicated expectation` = "Baseline", Gender = "-") %>% 
              bind_cols(baseline_data %>% 
                          distinct(Pigeon_ID) %>% 
                          summarise(`n pigeons observed` = length(Pigeon_ID)))) %>% 
  select(`Indicated expectation`, `n pigeons observed`, everything()) %>% 
  pander(split.cell = 20, split.table = Inf)

Indicated expectation	n pigeons observed	Gender	Estimated peck rate	Est.Error	Q2.5	Q97.5
Hungry	32	Woman	11.94	1.74	8.94	15.76
Hungry	12	Man	9.12	2.25	5.56	14.3
Satiated	25	Woman	7.88	1.34	5.61	10.85
Satiated	8	Man	7.34	2.25	4.05	12.78
Baseline	35	-	1.32	0.51	0.54	2.53

\(~\)

Build Figure 2e-h

\(~\)

Get posterior means and difference contrasts

# treatment model

peck_draws_treatment <-
  as_draws_df(peck_model_treatment) %>% 
  mutate(Hungry = exp(b_Bias_treatmentHungry),
         Satiated = exp(b_Bias_treatmentSatiated),
         diff_contrast = (Hungry - Satiated)) %>% 
  select(Hungry, Satiated, diff_contrast) %>%
  bind_cols(baseline_peck_predictions %>% select(Baseline_estimate)) %>% 
  pivot_longer(names_to = "Treatment", values_to = "Posterior_estimate", cols = 1:3) %>% 
  mutate(Predictor = "Allocated primer")


p6 <- 
    peck_draws_treatment %>% 
    filter(Treatment != "diff_contrast") %>% 
    ggplot(aes(x = Treatment, y = Posterior_estimate)) +
    stat_slab(aes(y = Baseline_estimate),
              linetype = 2, linewidth = 0.8, slab_fill = "white",
              colour = "black") +
    stat_halfeye(aes(fill = Treatment), .width = c(0.66, 0.95), alpha = 0.9,
                 point_interval = "median_qi", point_fill = "white",
                 shape = 21, point_size = 4, stroke = 1.5) +
    scale_fill_manual(values = met.brewer("Hiroshige", 2)) +
    coord_flip(ylim = c(0, 20)) +
    #geom_vline(xintercept = 0, linetype = 2) +
    #scale_y_continuous(breaks = c(, 0, 1)) +
    xlab("Allocated primer") +
    ylab("Estimated pecks per min") +
    theme_bw() + 
    theme(legend.position = "none",
          panel.grid.minor = element_blank(),
          text = element_text(size = 14))

p7 <-
  peck_draws_treatment %>% 
  filter(Treatment == "diff_contrast") %>% 
  ggplot(aes(y = Posterior_estimate)) + 
   stat_halfeye(aes(fill = Treatment), .width = c(0.66, 0.95), alpha = 0.9,
               point_interval = "median_qi", point_fill = "white",
               shape = 21, point_size = 4, stroke = 1.5, scale =0.5) +
  scale_fill_manual(values = met.brewer("Hiroshige")[4]) +
  coord_flip(ylim = c(-5, 12)) +
  geom_hline(yintercept = 0, linetype = 2, linewidth = 0.75) +
  #scale_y_continuous(breaks = c(, 0, 1)) +
  xlab(NULL) +
  ylab("Hungry - Satiated difference\ncontrast (pecks per min)") +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        text = element_text(size = 14))

# expectation model

peck_draws_expectation <-
  as_draws_df(peck_model_expectation) %>% 
  mutate(Hungry = exp(b_ExpectationHungry),
         Satiated = exp(b_ExpectationSatiated),
         diff_contrast = (Hungry - Satiated)) %>% 
  select(Hungry, Satiated, diff_contrast) %>% 
  bind_cols(baseline_peck_predictions %>% select(Baseline_estimate)) %>% 
  pivot_longer(names_to = "Treatment", values_to = "Posterior_estimate", cols = 1:3) %>% 
  mutate(Predictor = "Allocated primer")

p8 <-
  peck_draws_expectation %>% 
  filter(Treatment != "diff_contrast") %>% 
  ggplot(aes(x = Treatment, y = Posterior_estimate)) + 
  stat_slab(aes(y = Baseline_estimate),
              linetype = 2, linewidth = 0.8, slab_fill = "white",
              colour = "black") +
    stat_halfeye(aes(fill = Treatment), .width = c(0.66, 0.95), alpha = 0.9,
               point_interval = "median_qi", point_fill = "white",
               shape = 21, point_size = 4, stroke = 1.5) +
  scale_fill_manual(values = met.brewer("Hiroshige", 2)) +
  coord_flip(ylim = c(0, 20)) +
  #geom_vline(xintercept = 0, linetype = 2) +
  #scale_y_continuous(breaks = c(, 0, 1)) +
  xlab("Indicated expectation") +
  ylab("Estimated pecks per min") +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank(),
        text = element_text(size = 14))

p9 <-
  peck_draws_expectation %>% 
  filter(Treatment == "diff_contrast") %>% 
  ggplot(aes(y = Posterior_estimate)) + 
   stat_halfeye(aes(fill = Treatment), .width = c(0.66, 0.95), alpha = 0.9,
               point_interval = "median_qi", point_fill = "white",
               shape = 21, point_size = 4, stroke = 1.5, scale =0.5) +
  scale_fill_manual(values = met.brewer("Hiroshige")[4]) +
  coord_flip(ylim = c(-5, 12)) +
  geom_hline(yintercept = 0, linetype = 2, linewidth = 0.75) +
  #scale_y_continuous(breaks = c(, 0, 1)) +
  xlab(NULL) +
  ylab("Hungry - Satiated difference\ncontrast (pecks per min)") +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        text = element_text(size = 14))

Now for the gender models

# with gender

gender_treatment_draws_2 <-
  fitted(peck_model_treatment_gender, 
         newdata = new_data_gender, summary = F, re_formula = NA) %>% 
  as_tibble() %>% 
  rename(Hungry_Women = V1, Hungry_Men = V2, Satiated_Women = V3, Satiated_Men = V4) %>% 
  bind_cols(baseline_peck_predictions %>% select(Baseline_estimate)) %>% 
  pivot_longer(names_to = "Group", values_to = "Posterior_estimate", cols = 1:4) %>% 
  separate(sep = "_", col = Group, into = c("Treatment", "Gender"))
  
calculate_all_the_diffs_3 <-
  fitted(peck_model_treatment_gender, 
         newdata = new_data_gender, summary = F, re_formula = NA) %>% 
  as_tibble() %>% 
  rename(Hungry_Women = V1, Hungry_Men = V2, Satiated_Women = V3, Satiated_Men = V4) %>% 
  mutate(Women_h_s_diff = Hungry_Women - Satiated_Women,
         Men_h_s_diff = Hungry_Men - Satiated_Men,
         diff_diff = Women_h_s_diff - Men_h_s_diff) %>% 
  select(contains("diff")) %>%
  rename(`H-S (women)` = Women_h_s_diff,
         `H-S (men)` = Men_h_s_diff,
         `Interaction` = diff_diff) %>% 
  pivot_longer(cols = 1:3, names_to = "diff_contrast", values_to = "posterior_diff")

gp5 <- 
  gender_treatment_draws_2 %>% 
  ggplot(aes(x = Gender, y = Posterior_estimate)) + 
    stat_slab(alpha = 0.8, shape = 21, aes(fill = Treatment)) +
  stat_slab(aes(y = Baseline_estimate),
              linetype = 2, linewidth = 0.8, slab_fill = "white",
              colour = "black") +
  #stat_halfeye(.width = c(0.66, 0.95), alpha = 0.9,
   #            point_interval = "median_qi", point_fill = "white",
   #            shape = 21, point_size = 4, stroke = 1.5,
   #            fill = met.brewer("Hiroshige", 5)[2]) +
  scale_fill_manual(values = c(met.brewer("Hiroshige", 10)[4], met.brewer("Hiroshige", 10)[6])) +
  coord_flip()+#ylim = c(25, 55)) +
  #geom_vline(xintercept = 0, linetype = 2) +
  #scale_y_continuous(breaks = c(, 0, 1)) +
  labs(x = "Gender", y = "Estimated pecks per min", fill = "Allocated\nPrimer") +
  theme_bw() + 
    theme(legend.position = "bottom",
        panel.grid.minor = element_blank(),
        text = element_text(size = 12),
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 10))



gp6 <-
  calculate_all_the_diffs_3 %>% 
  ggplot(aes(x = diff_contrast, y = posterior_diff)) + 
    #stat_slab(alpha = 0.9, shape = 21, aes(fill = Gender)) +
  stat_halfeye(.width = c(0.66, 0.95), alpha = 0.9,
               point_interval = "median_qi", point_fill = "white",
               shape = 21, point_size = 3, stroke = 1.5,
               fill = met.brewer("Hiroshige", 10)[5]) +
  #scale_fill_manual(values = met.brewer("Hiroshige", 2)) +
  coord_flip()+#ylim = c(25, 55)) +
   geom_hline(yintercept = 0, linetype = 2, linewidth = 0.75) +
  #geom_vline(xintercept = 0, linetype = 2) +
  scale_y_continuous(limits = c(-10, 15)) +
  xlab("Difference contrast") +
  ylab("Pecks per min") +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank(),
        text = element_text(size = 12))


# expectation model

gender_expectation_draws_2 <-
  fitted(peck_model_expectation_gender, 
         newdata = new_data_gender_2, summary = F, re_formula = NA) %>% 
  as_tibble() %>% 
  rename(Hungry_Women = V1, Hungry_Men = V2, Satiated_Women = V3, Satiated_Men = V4) %>% 
  bind_cols(baseline_peck_predictions %>% select(Baseline_estimate)) %>% 
  pivot_longer(names_to = "Group", values_to = "Posterior_estimate", cols = 1:4) %>% 
  separate(sep = "_", col = Group, into = c("Treatment", "Gender"))
  
calculate_all_the_diffs_4 <-
  fitted(peck_model_expectation_gender, 
         newdata = new_data_gender_2, summary = F, re_formula = NA) %>% 
  as_tibble() %>% 
  rename(Hungry_Women = V1, Hungry_Men = V2, Satiated_Women = V3, Satiated_Men = V4) %>% 
  mutate(Women_h_s_diff = Hungry_Women - Satiated_Women,
         Men_h_s_diff = Hungry_Men - Satiated_Men,
         diff_diff = Women_h_s_diff - Men_h_s_diff) %>% 
  select(contains("diff")) %>%
  rename(`H-S (women)` = Women_h_s_diff,
         `H-S (men)` = Men_h_s_diff,
         `Interaction` = diff_diff) %>% 
  pivot_longer(cols = 1:3, names_to = "diff_contrast", values_to = "posterior_diff")

gp7 <- 
  gender_expectation_draws_2 %>% 
  ggplot(aes(x = Gender, y = Posterior_estimate)) + 
    stat_slab(alpha = 0.9, shape = 21, aes(fill = Treatment)) +
  stat_slab(aes(y = Baseline_estimate),
              linetype = 2, linewidth = 0.8, slab_fill = "white",
              colour = "black") +
  #stat_halfeye(.width = c(0.66, 0.95), alpha = 0.9,
   #            point_interval = "median_qi", point_fill = "white",
   #            shape = 21, point_size = 4, stroke = 1.5,
   #            fill = met.brewer("Hiroshige", 5)[2]) +
  scale_fill_manual(values = met.brewer("Hiroshige", 2)) +
  coord_flip()+#ylim = c(25, 55)) +
  #geom_vline(xintercept = 0, linetype = 2) +
  #scale_y_continuous(breaks = c(, 0, 1)) +
  labs(x = "Gender", y = "Estimated pecks per min", fill = "Indicated\nexpectation") +
  theme_bw() + 
  theme(legend.position = "bottom",
        panel.grid.minor = element_blank(),
        text = element_text(size = 12),
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 10))



gp8 <-
  calculate_all_the_diffs_4 %>% 
  ggplot(aes(x = diff_contrast, y = posterior_diff)) + 
    #stat_slab(alpha = 0.9, shape = 21, aes(fill = Gender)) +
  stat_halfeye(.width = c(0.66, 0.95), alpha = 0.9,
               point_interval = "median_qi", point_fill = "white",
               shape = 21, point_size = 3, stroke = 1.5,
               fill = met.brewer("Hiroshige", 10)[5]) +
  #scale_fill_manual(values = met.brewer("Hiroshige", 2)) +
  coord_flip()+#ylim = c(25, 55)) +
   geom_hline(yintercept = 0, linetype = 2, linewidth = 0.75) +
  #geom_vline(xintercept = 0, linetype = 2) +
  scale_y_continuous(limits = c(-10, 15)) +
  xlab("Difference contrast") +
  ylab("Pecks per min") +
  theme_bw() + 
  theme(legend.position = "none",
        panel.grid.minor = element_blank(),
        text = element_text(size = 12))

\(~\)

Find the difference between baseline and observer estimates

Table S5. The degree to which each group of observer’s overestimates feeding rate (number of ground pecks per minute)

baseline_peck_predictions %>% select(Baseline_estimate) %>%  bind_cols(
  
  as_draws_df(peck_model_treatment) %>% 
    mutate(Hungry = exp(b_Bias_treatmentHungry),
           Satiated = exp(b_Bias_treatmentSatiated)) %>% 
    select(Hungry, Satiated)) %>% 
  mutate(`Bias treatment Satiated / Baseline` = Satiated / Baseline_estimate,
         `Bias treatment Hungry / Baseline` = Hungry / Baseline_estimate) %>% 
  select(contains("Bias")) %>% 
  pivot_longer(cols = everything(), values_to = "estimate", names_to = "Stat") %>% 
  group_by(Stat) %>% 
  summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE), .cores = 4) %>% 
  select(-variable) %>% 
  
  bind_rows(
    
    baseline_peck_predictions %>% select(Baseline_estimate) %>%  bind_cols(
      
      as_draws_df(peck_model_expectation) %>% 
        mutate(Hungry = exp(b_ExpectationHungry),
               Satiated = exp(b_ExpectationSatiated)) %>% 
        select(Hungry, Satiated)) %>% 
      mutate(`Expectation Satiated / Baseline` = Satiated / Baseline_estimate,
             `Expectation Hungry / Baseline` = Hungry / Baseline_estimate) %>% 
      select(contains("Expectation")) %>% 
      pivot_longer(cols = everything(), values_to = "estimate", names_to = "Stat") %>% 
      group_by(Stat) %>% 
      summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE), .cores = 4) %>% 
      select(-variable)
  ) %>% 
  pander()

Stat	median	sd	2.5%	97.5%
Bias treatment Hungry / Baseline	7.883	4.111	3.728	19.04
Bias treatment Satiated / Baseline	6.713	3.562	3.149	16.48
Expectation Hungry / Baseline	8.67	4.554	4.108	20.94
Expectation Satiated / Baseline	5.722	3.065	2.688	14

Figure 2

Option 1

(p2 + p3) / (p4 + p5) / (p6 + p7) / (p8 + p9) +
   plot_annotation(tag_levels = 'a')

Option 2

(gp1 + gp2) / (gp3 + gp4) / (gp5 + gp6) / (gp7 + gp8) +
   plot_annotation(tag_levels = 'a')

Figure 2. Posterior mean estimates and difference contrasts for observer estimated group foraging percentage and individual feeding rates. The coloured area is the posterior distribution and the white point is the mean estimate with associated 67% and 95% credible intervals.

Session info

sessionInfo() %>% pander

R version 4.2.2 (2022-10-31 ucrt)

Platform: x86_64-w64-mingw32/x64 (64-bit)

locale: LC_COLLATE=English_United Kingdom.utf8, LC_CTYPE=English_United Kingdom.utf8, LC_MONETARY=English_United Kingdom.utf8, LC_NUMERIC=C and LC_TIME=English_United Kingdom.utf8

attached base packages: grid, stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: ggdag(v.0.2.8), png(v.0.1-8), pander(v.0.6.5), DT(v.0.28), patchwork(v.1.1.2), kableExtra(v.1.3.4), MetBrewer(v.0.2.0), tidybayes(v.3.0.4), loo(v.2.6.0), brms(v.2.19.0), Rcpp(v.1.0.10), lubridate(v.1.9.2), forcats(v.1.0.0), stringr(v.1.5.0), dplyr(v.1.1.2), purrr(v.1.0.1), readr(v.2.1.4), tidyr(v.1.3.0), tibble(v.3.2.1), ggplot2(v.3.4.2) and tidyverse(v.2.0.0)

loaded via a namespace (and not attached): backports(v.1.4.1), systemfonts(v.1.0.4), plyr(v.1.8.8), igraph(v.1.4.2), svUnit(v.1.0.6), crosstalk(v.1.2.0), rstantools(v.2.3.1), inline(v.0.3.19), digest(v.0.6.31), htmltools(v.0.5.5), viridis(v.0.6.3), fansi(v.1.0.4), magrittr(v.2.0.3), checkmate(v.2.2.0), tzdb(v.0.4.0), graphlayouts(v.1.0.0), RcppParallel(v.5.1.7), matrixStats(v.0.63.0), vroom(v.1.6.3), xts(v.0.13.1), svglite(v.2.1.1), timechange(v.0.2.0), prettyunits(v.1.1.1), colorspace(v.2.1-0), ggrepel(v.0.9.3), rvest(v.1.0.3), ggdist(v.3.3.0), xfun(v.0.39), callr(v.3.7.3), crayon(v.1.5.2), jsonlite(v.1.8.4), zoo(v.1.8-12), glue(v.1.6.2), polyclip(v.1.10-4), gtable(v.0.3.3), webshot(v.0.5.4), V8(v.4.3.0), distributional(v.0.3.2), pkgbuild(v.1.4.0), rstan(v.2.26.13), abind(v.1.4-5), scales(v.1.2.1), mvtnorm(v.1.1-3), miniUI(v.0.1.1.1), viridisLite(v.0.4.2), xtable(v.1.8-4), bit(v.4.0.5), stats4(v.4.2.2), StanHeaders(v.2.26.25), htmlwidgets(v.1.6.2), httr(v.1.4.6), threejs(v.0.3.3), arrayhelpers(v.1.1-0), posterior(v.1.4.1), ellipsis(v.0.3.2), pkgconfig(v.2.0.3), farver(v.2.1.1), sass(v.0.4.6), utf8(v.1.2.3), labeling(v.0.4.2), tidyselect(v.1.2.0), rlang(v.1.1.1), reshape2(v.1.4.4), later(v.1.3.1), munsell(v.0.5.0), dagitty(v.0.3-1), tools(v.4.2.2), cachem(v.1.0.8), cli(v.3.5.0), generics(v.0.1.3), evaluate(v.0.21), fastmap(v.1.1.1), yaml(v.2.3.7), processx(v.3.8.1), knitr(v.1.42), bit64(v.4.0.5), tidygraph(v.1.2.3), ggraph(v.2.1.0), nlme(v.3.1-160), mime(v.0.12), xml2(v.1.3.4), compiler(v.4.2.2), bayesplot(v.1.10.0), shinythemes(v.1.2.0), rstudioapi(v.0.14), curl(v.5.0.0), tweenr(v.2.0.2), bslib(v.0.4.2), stringi(v.1.7.12), highr(v.0.10), ps(v.1.7.5), Brobdingnag(v.1.2-9), lattice(v.0.20-45), Matrix(v.1.5-1), markdown(v.1.7), shinyjs(v.2.1.0), tensorA(v.0.36.2), vctrs(v.0.6.2), pillar(v.1.9.0), lifecycle(v.1.0.3), jquerylib(v.0.1.4), bridgesampling(v.1.1-2), httpuv(v.1.6.10), R6(v.2.5.1), promises(v.1.2.0.1), gridExtra(v.2.3), codetools(v.0.2-18), boot(v.1.3-28), colourpicker(v.1.2.0), MASS(v.7.3-58.1), gtools(v.3.9.4), withr(v.2.5.0), shinystan(v.2.6.0), parallel(v.4.2.2), hms(v.1.1.3), coda(v.0.19-4), rmarkdown(v.2.21), ggforce(v.0.4.1), shiny(v.1.7.4), base64enc(v.0.1-3) and dygraphs(v.1.1.1.6)